/** @addtogroup dft

 *  @{

 */

/*

  Copyright (C) 2016 D Levin (https://www.kfrlib.com)

  This file is part of KFR

  KFR is free software: you can redistribute it and/or modify

  it under the terms of the GNU General Public License as published by

  the Free Software Foundation, either version 3 of the License, or

  (at your option) any later version.

  KFR is distributed in the hope that it will be useful,

  but WITHOUT ANY WARRANTY; without even the implied warranty of

  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License

  along with KFR.

  If GPL is not suitable for your project, you must purchase a commercial license to use KFR.

  Buying a commercial license is mandatory as soon as you develop commercial activities without

  disclosing the source code of your own applications.

  See https://www.kfrlib.com for details.

 */

#pragma once

#include "../base/complex.hpp"

#include "../base/constants.hpp"

#include "../base/digitreverse.hpp"

#include "../base/read_write.hpp"

#include "../base/sin_cos.hpp"

#include "../base/small_buffer.hpp"

#include "../base/univector.hpp"

#include "../base/vec.hpp"

#include "../base/memory.hpp"

#include "../data/sincos.hpp"

CMT_PRAGMA_MSVC(warning(push))

CMT_PRAGMA_MSVC(warning(disable : 4127))

namespace kfr

{

namespace internal

{

template <typename T, size_t N, KFR_ENABLE_IF(N >= 2)>

CMT_INLINE vec<T, N> cmul_impl(const vec<T, N>& x, const vec<T, N>& y)

{

    return subadd(x * dupeven(y), swap<2>(x) * dupodd(y));

}

template <typename T, size_t N, KFR_ENABLE_IF(N > 2)>

CMT_INLINE vec<T, N> cmul_impl(const vec<T, N>& x, const vec<T, 2>& y)

{

    vec<T, N> yy = resize<N>(y);

    return cmul_impl(x, yy);

}

template <typename T, size_t N, KFR_ENABLE_IF(N > 2)>

CMT_INLINE vec<T, N> cmul_impl(const vec<T, 2>& x, const vec<T, N>& y)

{

    vec<T, N> xx = resize<N>(x);

    return cmul_impl(xx, y);

}

/// Complex Multiplication

template <typename T, size_t N1, size_t N2>

CMT_INLINE vec<T, const_max(N1, N2)> cmul(const vec<T, N1>& x, const vec<T, N2>& y)

{

    return internal::cmul_impl(x, y);

}

template <typename T, size_t N, KFR_ENABLE_IF(N >= 2)>

CMT_INLINE vec<T, N> cmul_conj(const vec<T, N>& x, const vec<T, N>& y)

{

    return swap<2>(subadd(swap<2>(x) * dupeven(y), x * dupodd(y)));

}

template <typename T, size_t N, KFR_ENABLE_IF(N >= 2)>

CMT_INLINE vec<T, N> cmul_2conj(const vec<T, N>& in0, const vec<T, N>& in1, const vec<T, N>& tw)

{

    return (in0 + in1) * dupeven(tw) + swap<2>(cnegimag(in0 - in1)) * dupodd(tw);

}

template <typename T, size_t N, KFR_ENABLE_IF(N >= 2)>

CMT_INLINE void cmul_2conj(vec<T, N>& out0, vec<T, N>& out1, const vec<T, 2>& in0, const vec<T, 2>& in1,

                           const vec<T, N>& tw)

{

    const vec<T, N> twr   = dupeven(tw);

    const vec<T, N> twi   = dupodd(tw);

    const vec<T, 2> sum   = (in0 + in1);

    const vec<T, 2> dif   = swap<2>(negodd(in0 - in1));

    const vec<T, N> sumtw = resize<N>(sum) * twr;

    const vec<T, N> diftw = resize<N>(dif) * twi;

    out0 += sumtw + diftw;

    out1 += sumtw - diftw;

}

template <typename T, size_t N, KFR_ENABLE_IF(N > 2)>

CMT_INLINE vec<T, N> cmul_conj(const vec<T, N>& x, const vec<T, 2>& y)

{

    vec<T, N> yy = resize<N>(y);

    return cmul_conj(x, yy);

}

template <typename T, size_t N, KFR_ENABLE_IF(N > 2)>

CMT_INLINE vec<T, N> cmul_conj(const vec<T, 2>& x, const vec<T, N>& y)

{

    vec<T, N> xx = resize<N>(x);

    return cmul_conj(xx, y);

}

template <typename T, size_t N>

using cvec = vec<T, N * 2>;

template <size_t N, bool A = false, typename T>

CMT_INLINE cvec<T, N> cread(const complex<T>* src)

{

    return cvec<T, N>(ptr_cast<T>(src), cbool_t<A>());

}

template <size_t N, bool A = false, typename T>

CMT_INLINE void cwrite(complex<T>* dest, const cvec<T, N>& value)

{

    value.write(ptr_cast<T>(dest));

}

template <size_t count, size_t N, size_t stride, bool A, typename T, size_t... indices>

CMT_INLINE cvec<T, count * N> cread_group_impl(const complex<T>* src, csizes_t<indices...>)

{

    return concat(read<N * 2, A>(ptr_cast<T>(src + stride * indices))...);

}

template <size_t count, size_t N, size_t stride, bool A, typename T, size_t... indices>

CMT_INLINE void cwrite_group_impl(complex<T>* dest, const cvec<T, count * N>& value, csizes_t<indices...>)

{

    swallow{ (write<A>(ptr_cast<T>(dest + stride * indices), slice<indices * N * 2, N * 2>(value)), 0)... };

}

template <size_t count, size_t N, bool A, typename T, size_t... indices>

CMT_INLINE cvec<T, count * N> cread_group_impl(const complex<T>* src, size_t stride, csizes_t<indices...>)

{

    return concat(read<N * 2, A>(ptr_cast<T>(src + stride * indices))...);

}

template <size_t count, size_t N, bool A, typename T, size_t... indices>

CMT_INLINE void cwrite_group_impl(complex<T>* dest, size_t stride, const cvec<T, count * N>& value,

                                  csizes_t<indices...>)

{

    swallow{ (write<A>(ptr_cast<T>(dest + stride * indices), slice<indices * N * 2, N * 2>(value)), 0)... };

}

template <size_t count, size_t N, size_t stride, bool A = false, typename T>

CMT_INLINE cvec<T, count * N> cread_group(const complex<T>* src)

{

    return cread_group_impl<count, N, stride, A>(src, csizeseq_t<count>());

}

template <size_t count, size_t N, size_t stride, bool A = false, typename T>

CMT_INLINE void cwrite_group(complex<T>* dest, const cvec<T, count * N>& value)

{

    return cwrite_group_impl<count, N, stride, A>(dest, value, csizeseq_t<count>());

}

template <size_t count, size_t N, bool A = false, typename T>